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Abstract 

o 

(N 

' We study the coupling of a tetraquark system to an exchanged meson-meson channel, using a pure 

gluonic theory based four-quark potential matrix model which is known to fit well a large number of 
data points for lattice simulations of different geometries of a four-quark system. We find that if this 
minimal-area-based potential matrix replaces the earlier used simple Gaussian form for the gluon field 
overlap factor / in its off-diagonal terms, the resulting T-matrix and phase shifts develop an angle 
dependence whose partial wave analysis reveals D wave and higher angular momentum components 
in it. In addition to the obvious implications of this result for the meson- meson scattering, this new 
O |' feature indicates the possibility of orbital excitations influencing properties of meson-meson molecules 

through a polarization potential. We have used a formalism of the resonating group method, treated 
kinetic energy and overlap matrices on model of the potential matrix, but decoupled the resulting 
complicated integral equations through the Born approximation. In this exploratory study we have 
used a quadratic confinement and not included the spin-dependence; we also used the approximation 
of equal constituent quark masses. 

in 

> , 

^ ■ 1 Introduction 

, Models of hadronic physics are to be compared with experimental results as well as with our understand- 

ing of QCD for large momentum transfers and for low momentum transfers where Feynman diagrams are 
not useful. In this way, models may help us in improving our knowledge of QCD for low and interme- 
diate energies of interest to hadronic spectroscopy and eventually nuclear physics. One way to get this 
understanding is to note some features present in the perturbative QCD, lattice gauge theory or models 
of atomic and nuclear physics and check if these features can be used in the low and intermediate energy 
hadronic physics. 

One such feature is an approach based on pair- wise interaction for an interacting multiparticle system 
(composed of more than two or three particles). This has been successful in atomic and many-nucleon 
systems; the corresponding two-body interaction being described by Coulombic and Yukawa potential, 
for example. The question is if the explicit presence of Non-Abelian gluon field can also be replaced 
by two body interquark potentials. The simplest way to use such a model is to try a sum of two-body 
potentials or interactions, the usual approach of atomic and nuclear physics. For comparison, it can be 
noted that the lowest order perturbative Feynman diagrams amplitudes are of this form, and a simple 
extension of this diagrammatic approach to multiquarks also has this pattern; see ref. [T] , and the later 
ones in the same approach, where the one gluon exchange potential, though, is replaced by Coulombic- 
plus-linear-plus-hypefine. If the numerically calculated energies of the four-quark systems on a lattice, 
in the static quark limit, are compared with a model that use only a sum of two-quark potentials, the 
model give a gross overestimate of the (magnitude of) four-quark binding energies; see fig. 4 of the same 
ref. [5]. A gluon field overlap factor / was introduced [3J essentially as a solution to this discrepancy. 
This factor multiplied only off-diagonal elements of the overlap, kinetic and potential energy matrices 
of the otherwise pairwise-sum-based Hamiltonian in the three-dimensional basis of the model system of 
four valence quarks/antiquarks and the purely gluonic field between them. Initially [3J the geometrical 
dependence of / on the quark positions was chosen purely based on computational convenience and had 
no known comparison with any QCD simulations. But when its different forms were compared [3] with 
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two-colour lattice numerical simulations in the pure gluonic theory, a factor of exp(— kS m i n ) had to be 
included in off-diagonal terms of the every version of the potential and overlap matrices; S m in is the 
minimal spatial area bounded by four external lines joining the four quarks and k is geometrically a 
constant. Only this way, a version of the model was eventually able to fit well "100 pieces of data — 
the ground and first excited states of configurations from six (kinds of) different four-quark geometries 
(including squares, rectangles, quadrilaterals and tetrahedra) calculated on a 16 3 x 32 lattice — with only 
four independent parameters to describe the interactions connecting the basis states" . It is to be noted 
that this exponential dependence on the spatial area in the model can be possibly traced back to the 
related use of the space-time area in more familiar models of Wilson loops studying time evolutions 
of a quark-antiquark pair. The connection was first suggested by a Wilson loop matrix in a strong 
coupling expansion scheme (see figs. 4 and 5 of ref. [5]) and appeared in above mentioned model [4] of 
the numerically evaluated Wilson loop matrix of the SU(2) C simulations. A full SU(3) C simulation [6] 
performed a bit later again showed the need for / model, though detailed geometrical dependence of / 
could not be studied in this 3 colour lattice gauge study. But a later numerical lattice study [7] by a 
reputed group of the full 2x2 matrix of the Wilson loops correlators and of " the interaction energy of the 
confining strings in the static rectangular tetraquark system in SU(3) gluodynamics" was well modeled 
again by a surface model, namely their soap film model that also incorporates (a flip to) the multi Y type 
linear potential emerging from recent numerical simulations [8] . The basis state overlap g [7] in the soap 
film model has a role similar to the gluon field overlap factor / of ref. [4] ; both / and g appear only in 
the off-diagonal terms of the respective matrices (N and T) of overlaps of the basis states. Continuing 
on the spatial and space-time areas, it can be pointed out that both kind of areas appear in eq.(I3) of 
the ref. [7] and are related to Wilson loops as earlier eq.(12) there indicates. 

Actually, above are models of the matrices of pure gluonic theory Wilson loops. The diagonal terms 
in these matrices are time evolutions of a tetraquark clustering and off-diagonal terms [7] are for time 
evolution that start from one tetraquark clustering (or topology) and end at another one. The numerical 
evaluation of the off-diagonal Wilson loops has been perhaps done in refs. [1] (and previous ones by the 
same group) and [7] only. But the diagonal Wilson loops have been studied by many other groups, most 
familiar being the studies reported in ref. [8] and the previous works by the same group. For one set of 
quark configurations, the diagonal studies are limited to only one Wilson loop. In a sense, this means 
limitation to only one state of the gluonic field as well, namely the one with the least energy to which 
system flips; if there are other states mentioned in the literature, these arc cither for comparison purpose 
{from which the system flips) or the excited state "contaminations". But a general study should actually 
incorporate a variety of basis states and thus include off-diagonal Wilson loops as well. The spatial area 
we are working on appears only in the off-diagonal Wilson loops and thus our work is not to be confused 
with the usual study of the diagonal Wilson loops effects. It is to be admitted that works like ref. [8] 
have indicated improvements in both evaluations and models of the diagonal Wilson loops and we have 
not included these improvements in our model of the diagonal term. But this is not a serious flaw, as 
a dynamical study [9] using these improved diagonal models mentions in its conclusions and comments 
that the " dynamics of (tetraquark) binding is dominated by the simple flip-flop term" , meaning that the 
essentially new connected string (butterfly) term introduced through the work of ref. [8] is dynamically 
" not rewarding" . It is to be noted that our diagonal terms include both terms whose minimum is the 
flip-flop term. 

This advocates the exp(— kS) form of / for static two quarks and two antiquarks. For a comparison 
with actual (hadron) experiments, we have to incorporate quark motion as well, possibly through using 
quark wave functions. The resulting four-body Schrodinger equation can be solved, as in ref. [TO] , 
variationally for the ground state of the system and the effective meson- meson potentials. Alternatively, 
the Hamiltonian emerging from the q 2 q 2 model has been diagonalized in the simple harmonic oscillator 
basis [llj . or was sandwiched between the external meson wave functions to give a transition amplitude 
of the Born diagrams [TJ 112) that is related to meson-meson phase shifts. We have used a formalism 
(resonating group method [13j ) that was, for the q 2 q 2 system, originally [3] used in a way that allowed 
finding inter-meson dependence with the dependence on a quark-antiquark distance being pre-fixed. 
The formalism allows using the best available knowledge of a meson wave function, though a simple 
Gaussian form for the wave functions and correspondingly a quadratic quark-antiquark potential was 
used for computational convenience. We have used this same formalism that can be generalized to using 
realistic meson wave functions and finding the inter-cluster dependence. But because of the additional 
computational problems due to a totally non-separable exponential of (a negative constant times) area 
in /, presently we had to also pre-specify a plane- wave form of the intcr-cluster dependence along with 



2 



still using Gaussian wave functions; effectively this means using a Born approximation as well. We have 
pointed out, though, that these wave function approximations are better than what their naive impression 
conveys: The Gaussian dependence on the quark- antiquark distance has the potential of resembling the 
realistic meson wave functions through an adjustment of its parameters [14] . And the plane wave form of 
the inter-cluster dependence is justified through a feeble inter-cluster (meson-meson) interaction noted in 
previous works [101 |3l [1] ; the meson-meson phase shifts resulting from this work are also much less than a 
radian. Only by using the Born- approximation, the resulting coupled integral equations for inter-cluster 
wave functions could be decoupled in this work. This decoupling allowed us to numerically calculate 
the off-diagonal elements as nine-dimensional integrals for the components of the eventual four position 
3-vectors; only the overall center-of-mass dependence could be analytically dealt with in a trivial manner. 
Before this numerical integration, for the kinetic energy terms we had to differentiate the area in /. The 
form of area used in the detailed form of the Q 2 Q 2 model, that we take from ref. [4], has square roots 
of the functions of our position variables. Thus a differentiation of this area form yields in denominators 
combinations of position variables that become zero somewhere in the ranges of integrations to be later 
done. The numerical evaluations of the resulting nine-dimensional improper integrals is expected to be too 
demanding, as our initial explorations indicated. Thus, for the to-be-differentiated right y 7 / part of the 
some kinetic energy terms we replaced the area by an approximated quadratic form whose differentiation 
does not result in negative powers of the position variables. 

We also find that the use of the / factor in the new form reduces the long range meson-meson 
interaction and thus as usual solves the well known Van der Walls force problem 115/91 with the otherwise 
naive sum of one gluon exchange pair-wise interaction. It has been said [10] that dynamically it is not a 
serious problem because of the quark-antiquark pair creation and because of the wave function damping 
of the large distance configurations. Though ref. [6 partly incorporated both the quark-antiquark pair 
creation and the meson wave functions and still showed a need for the / factor, through the present work 
we want to point out that the dynamical role of the / factor in meson-meson interactions is not limited to 
solving Van der Walls force problem or pointing out [16] otherwise over-binding in certain meson-meson 
systems. The / factor points certain features (non-separability of /) of QCD that are 1) indicated by 
lattice simulations and 2) can be compared with actual experiments. 

There have been recent hadron- level studies [17) [9] using the above mentioned improved models of 
the diagonal Wilson loops. But the quark level limitations of the models [H][8] mean similar limitations 
for the hadron-level results: [17] [5] study the properties (like binding energy and direct potential of 
the ground state itself (or in isolation) , whereas we aim to study the dynamical coupling of a tetraquark 
state to other basis state(s) of the tetraquark system — essentially to the other clustering or the exchanged 
channel. Thus, as we say in the abstract, in addition to doing the phase-shift calculations for a meson- 
meson scattering, we study a coupling that can affect even the ground state itself through a second 
order perturbation theory effect named polarization potential [19] , That is, after including the quark 
mass differences, a meson-meson state may not be degenerate with an exchanged channel and thus the 
coupling between this state and the exchanged intermediate one (a hadron loop) may help resolve the 
underlying structure of a possible meson-meson state. Such a state may be a meson-meson molecule that 
can be formed by the ground state. This also makes a study of the dynamical coupling of a meson-meson 
channel to exchanged one worth pursuing. 

In section 2 we have written the total state vector of the q 2 q 2 system as in RGM, along with introducing 
the Hamiltonian H of the system without the / factor and then modifying H through the /. In section 
3 different position dependent forms of / have been described, including the approximate forms that we 
had to use. In section 4 we have solved the integral equations for a meson meson molecule in the absence 
of spin degrees of freedom and with all equal quark masses. This section ends with a prescription to 
find the phase shifts. In the last section we have presented the numerical values of the phase shifts for 
different forms of /, for different values of free parameter kt and for different values of angle between 
Pi and P 2 . 

2 The Q 2 Q 2 Hamiltonian and the wave-function 

Using adiabatic approximation we can write the total state vector of a system containing two quarks two 
antiquarks and the gluonic field between them as a sum of the product of the quarks (Q or Q) position 
dependent function \E' g (ri, r2, rg, rj) and the gluonic field state |fc) fl . \k) g is defined as a state which 
approaches \k) c in the weak coupling limit , with |l) c = |ligl2i)o |2) c = |liil23)c and |3) c = |3i23g4) c . 
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In lattice simulations of the corresponding (gluonic) Wilson loops it is found that the lowest eigenvalue 
of the Wilson matrix, that is energy of the lowest state, is always the same for both 2x2 and 3x3 
matrices provided that |l) fl or \2) g has the lowest energy [4|. The later calculations [7] of the tetraquark 
system were also done with a two level approximation. Taking advantage of these observations, we have 
included in our expansion only two basis states. As in resonating group method, ^(ri, r 2 , r^, r^) or 
$ 9 (R C . Rfc, y fe , Zfe) is written as product of known dependence on R c ,y fe ,Zfc and unknown dependence 
on Rfc. i.e. 4 , 3 (ri, r 2 , rg, rj) = ^ c (R c )xfe(Rfc)V'fe(yfei z fc)- Here R c is the center of mass coordinate of 
the whole system, Ri is the vector joining the center of mass of the clusters (1,3) and (2,4), y 1 is the 
position vector of quark 1 with respect to 3 within the cluster (1, 3) and Zi is the position vector of quark 
2 with respect to 4 within the cluster (2,4). The same applies to R2, y 2 and z 2 for the clusters (1,4) and 
(2,3). Similarly we can define R3,y 3 and Z3 for the clusters (1,2) and (3,4). Or we can write them in 
terms of position vector of the four particles (quarks or antiquarks) as follow 

R4 = ^(ri + rj- r 2 - r ¥ ) ) y 1 =r 1 -r s and z x = r 2 - r ¥ , (1) 

R-2 = ^(ri + r ¥ - r 2 - rg-) , y 2 = ri - r ¥ and z 2 = r 2 - rg- (2) 

and 

R 3 = q ( r l + r 2 - r 3 ~ r i) .y3= r l- r 2 and z 3 = - Fj. (3) 

Thus meson meson state vector in the restricted gluonic basis is written as 

2 

|*(n,r 2 ,r3,r4; 5 )) = ^ |fc) ff * c (R c )xfc(Rfc)4(y fc )Cfe(zfe)- (4) 
fe=l 

2 2 

Here Ck(Yk) — — 5" ex P(^^) an d Cfc( z fc) = — 3- ex p(x7r)- These Gaussian forms of meson wave 

(27rd 2 )4 4a (27rd 2 )4 4a 

functions are, strictly speaking, the wave functions of a quadratic confining potential. But, as pointed 
out in text below fig. 1 of ref. |14) . the overlap of a Gaussian wave function and the eigenfunction of the 
realistic linear plus colombic potential can be made as close as 99.4% by properly adjusting its parameter 
d. A realistic value of d mimicking a realistic meson wave function depends on the chosen scattering 
mesons and thus is postponed to our future work |20j . Presently, to explore the qualitative implications 
of the geometric features of the gluonic overlap factor /, we have used in £fc(yj.) and Cfc( z fc) a value 
d = 0.558 fm defined by the relation d 2 = y/3R%/2 [TB], with R c = 0.6 fm being the r.m.s. charge radius 
of the qqq system whose wave function is derived by using the same quadratic confining potential. 
As for the Hamiltonian, for f—1 the total Hamiltonian H of our 4-particle system is taken as [21] 

4 pi 4 

6 = E h + 2^1 + £ u ( r *i) F * ,F J- (5) 

i=l ' i<j 

Our same constituent quark mass value m = 0.3GeV for all quarks and antiquarks is one used in refs. [3], 
and our kinetic energy operator is similarly non-relativistic; it is included in our aims to compare with 
this work and isolate the effects only due to a different expression for the /. In above each of Fj has 8 
components F- — A /2, I = 1,2, 3, 8 and Fj* = A'*/2, A' are Gell-Mann matrices operating on the i-th 
particle; I is shown as a superscript only to avoid any possible confusion with subscript i which labels a 
particle. 

For the pairwise qq potential, we have used a quadratic confinement 

Vij = Cr 2 %] + C with i, j = 1, 2, 3, 4. (6) 

for exploratory study. While we have neglected short range coulomb like interactions as well as spin- 
dependent terms. Along with that non relativistic limit has also been taken. The model used by Vijande 
[9] is also restricted to these limits. As for the within-a-cluster dependence of the wave function, this 
use of the quadratic potential in place of the realistic Coulumbic plus linear may change the full wave 
function. In the within-a-cluster, this change of wave function is found to result in a change of an overlap 
integral from 100% to 99.4% only provided the parameter d of the wave function is adjusted. Although the 
expression for P c written immediately after eq. (1261) suggests a way to connect the additional parameter of 
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the full wave function with that of a cluster, it is difficult to make a similar overlap test for the full wave 
function. But there is no a priori reason to deny that at least the qualitative features (like a new kind of 
angle dependence mentioned in the results part) we want to point out using this quadratic confinement 
would survive in a more realistic calculation; a similar exploration of the q 2 q 2 system properties was 
first done |10) using a quadratic confinement and later extended and improved calculation [21] with 
more realistic pair-wise interaction reinforced the KK results obtained initially through the quadratic 
confinement. It seems that a proper adjustment the parameters of the quadratic (or SHO) model can 
reasonably simulate a qq or even a q 2 q 2 system. In our case, this adjustment of the parameters can be 
done once a choice of actual scattering mesons is made in a formalism [20] that incorporates spin and 
flavour degrees of freedom. But, as shown in fig. 2(b) of ref.|10]. properties of the q 2 q 2 system may not 
be very sensitive to the actual values of the parameters and we expect our presently chosen values of the 
parameters to well indicate the essential features resulting from the non-separable form of the gluonic 
field overlap factor /. 

For the central simple harmonic oscillator potential of eq.(j6|), the above mentioned size d in the 
eigenfunctions £,k(yk) an d Cfc( z fc) is related to the quadratic coefficient C which thus is given a value of 
-0.0097GeV 3 . 

As in a resonating group calculation, we take only variations in the Xfc factor of the total state vector of 
the system written in eq.(j4|). Setting the coefficients of linearly independent arbitrary variations <5xfc(R.fc) 
as zero and integrating out R c , (Sip \ H — E c \ ip) = from eq.j4} gives 

I d 3 y k d 3 z k Uy k )U^) g (k \H-E C \ l) B xi(R*)bbiM*i) = 0, (7) 

for each of the k values (1 and 2). According to the (2 dimensional basis) model I a of ref. [I], the 
normalization, potential energy and kinetic energy matrices in the corresponding gluonic basis are 

W = (i/V). o 

y = ( TT ( W 13 + V 2l) |/(«ia + «S3 :_- V 13 - V 24 - v 14 - v 2s) 

V lf{vi2+m- V 13- V 24~ V 14~ V 2s) "T Ki + V 23 ) 

and 

g (k\ K \ o s =wA(E-3)^(/)w- ( 10 ) 

i=l 

This is the modification, through the / factor, to the Hamiltonian as much as we need it for the integral 
equations below in section 4 (that is only the modified matrix elements). 



3 Different forms of / 

Ref. [1] supports through a comparison with numerical lattice simulations a form of / that was earlier [15] 
suggested through a quark-string model extracted from the strong coupling lattice Hamiltonian gauge 
theory. This is 

/ = cxp(-b s k f S), (II) 

S being the area of minimal surface bounded by external lines joining the position of the two quarks and 
two antiquarks, and b s = 0.18GeV 2 is the standard string tension |22 [ ll8 j . fc/ is a dimensionless parameter 
whose value of 0.57 was decided in ref. [4] by a fit of the simplest two-state area-based model (termed 
model Ia) to the numerical results for a selection of Q 2 Q 2 geometries. It is shown there [1] that the 
parameters, including kf, extracted at this SU(2) C lattice simulation with (3 — 2.4 can be used directly 
in, for example, a resonating group calculation of a four quark model as the continuum limit is achieved 
for this value of /3. 

The simulations reported in ref. [1] were done in the 2-colour approximation. But, for calculating the 
dynamical effects, we use actual SU(3) colour matrix elements of ref. [3]. The only information we take 
from the computer simulations of ref. [1] is value of kf. This describes a geometrical property of the 
gluonic field (its spatial rate of decrease to zero) and it may be the case that the geometrical properties 
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of the gluonic field are not much different for different number of colours, as suggested for example by 
successes of the geometrical flux tube model. Situation is more clear, though, for the mass spectra and the 
string tension generated by the gluonic field: ref. [23] compare these quantities for SU(2) C , SU(3) C and 
SU(4) C gauge theories in 2+1 dimensions and find that the ratio of masses are, to a first approximation, 
are independent of the number of colours. Their preliminary calculations in 3+1 dimensions indicate 
a similar trend. Directly for the parameter kf, appearing in the overlap factor / studied in this work, 
the similar conclusion can be drawn from a comparison of the mentioned lattice calculations [6] on the 
interaction energy of the two heavy-light Q 2 q 2 mesons in the realistic 5£/(3) c gauge theory with ref. 
[4] that uses SU(2) C . For interpreting the results in terms of the potential for the corresponding single 
heavy-light meson (Qq), a Gaussian form 

f = eM-b s k f J2 r i j ) (12) 

i<j 

of the gluonic filed overlap factor / is used in this ref. [5] for numerical convenience and not the minimal 
area form. But for a particular geometry, the two exponents (the minimal area and the sum of squared 
distances) in these two forms of / are related and thus for a particular geometry a comparison of the 
parameter kf multiplying area and corresponding (different!) kf multiplying sum of squares in eq. (14) 
of ref. [6] is possible. We note that, after correcting for a ratio of 8 between the sum of distance squares 
(including two diagonals) and the area for the square geometry, the colour-number-generated relative 
difference for this geometry is just 5%: the coefficient is 0.075 x 8 = 0.6 multiplying sum of squared 
distances and 0.57 multiplying the minimal area. But, as the precise form of / is still under development 
(the latest work [7] has covered only a very limited selection of the positions of tetraquark constituents) 
and the expression for the area in its exponent needs improvement, it is not sure precisely what value 
of the kf best simulates QCD and we have mainly worked with an approximate value of 0.5 that is also 
mentioned in ref. |4| and is numerically easier to deal with. (It is to be noted that the soap film model of 
ref. [7] does not treat kf as a variational parameter. If that is interpreted as fixing kf at 1, this prescription 
might have been successful due to their selection of quark configurations being limited to planar ones; a 
work [24] by UKQCD that is limited to planar geometries also favors a value closer to 1. But their more 
general work [4] resulted in a value of kf near 0.5.) 

For the area as well, ref. [4] used an approximation: A good model of area of the minimal surface 
could be that given in ref. [25] as 

5 = du dvKur-fi + (1 - u)rj 2 ) x (vr^ + (1 -v)^)] (13) 
Jo Jo 

(Work is in progress [26) . to judge the surface used in this model, and its area, from the point of 
differential geometry and there are indications that this is quite close to a minimal surface.) But the 
simulations reported in ref. [4] were carried out for the 5 in eq. (|llj) being "the average of the sum of the 
four triangular areas defined by the positions of the four quarks" . Although for the tetrahedral geometry 
the 5 used in ref. [4] is as much as 26 percent larger than the corresponding minimal-like area of ref. 
|25j . it can be expected that their fitted value of kf is reduced to partially compensate this over estimate 
of the 5 area. Anyway, as we are calculating the dynamical effects of the model of ref. [4] , we have used 
the form of 5 that is used in this work. 

The area 5 of ref. [4] becomes (with a slight renaming) 



5 = -[5(134) + 5(234) + 5(123) + 5(124)], (14) 

where S(ijk) is the area of the triangle joining the vertices of the positions of the quarks labled as i,j 
and k. In the notation of eqs.([l]|3]) this becomes 5(134) = ^\y 1 x z 3 | = ||(R.2 + R3) x (Ri — Ra)!) 
5(234) = ±|z 2 x z 3 | = ±|(R 3 - Ri) x (R, - R 2 )|, 5(123) = ±|y 3 x z 2 | = ^(H* + R 2 ) x (R 3 - 
and 5(124) = i|y 3 X ZjJ = ||(Ri + R 2 ) X (R 3 — R 2 )|. Written in terms of the rectangular components 
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(xi,yi,zi) of Ri, {x 2 ,y2,z 2 ) of R 2 and (x 3 ,y 3 ,z 3 ) of R 3 , we have 

5(134) = ~ | (x 2 2 + 2/2 2 + z 2 2 + 2 (x 2 x 3 + y 2 y 3 + z 2 z 3 ) + x 3 2 + y 2 + z 3 2 ) 
(xi 2 + y 2 + z 2 - 2 (X1X2 + 2/12/2 + ziz 2 ) + x 2 2 + 2/2 2 + z 2 2 ^j 
- [xix 2 + 2/12/2 + z 1 z 2 - (x 2 2 + 2/2 2 + z 2 2 ) + x x x 3 + y x y 3 + z 1 z 3 - (x 2 x 3 + y 2 y 3 + z 2 z 3 ) ) j . (15) 

Explicit rectangular expressions for 5(234), 5(123) and 5(124) are similar. 

This form of 5 has square roots. For the K.E. part of the Hamiltonian matrix (see eq.fTOjl ), we have 
to differentiate an exponential of this square root. After differentiating we can have negative powers of 5 
and when they will be integrated in the latter stages can have singularities in the integrands resulting in 
computationally too demanding nine-dimensional improper integrals. Thus we have availed to ourselves 
an approximated 5, named S a , which is a sum of different quadratic combinations of quarks positions. We 
chose S a by minimizing J d 3 Rid 3 R,2d 3 R 3 (S — S a ) 2 with respect to the coefficients of the quadratic position 
combinations; that is, these coefficients are treated as variational parameters. The first (successful) form 
which we tried for S a was 

S a =a (xi 2 + 2/i 2 + zi 2 ) + b (x 2 2 + 2/ 2 2 + ^ 2 2 ) + C (a; 3 2 + y 3 2 + z 3 2 ) + d Xl x 2 + 

eyiV2 + f z\Z2 + gx 2 x 3 + hy 2 y 3 + iz 2 z 3 + jx x x 3 + ky%y 3 + lz\z 3 . (16) 

This contained 12 variational parameters a,b,c,...,l. Minimization gave values (reported with accuracy 4 
though in the computer program accuracy 16 was used)as 

a = 0.4065, 6 = 0.4050, c= 0.3931,j= -0.0002, 1 = -0.0002. 

In the reported accuracy other parameters are zero. Here limits of integration were from -15 to 15 in 
GeV -1 . We also tried S a as 

3 3 3 

ij^iVj ~l~ Ci,j%i z j + dijyiZj 

i=l *,j=l i<j,j=2 

and 

3 3 3 

^ i/.-'T + I",!/; + riiz 2 ) + ^ ( b i,3 x iVj + £i,j x i z j + <k,jViZj) + ^ ( e i,o x i x i + fijViVi + 9i>jZiZj), 

i=l i,j=l i<j,j=2 

with variational parameters being 39 and 45 respectively. Both the latter forms gave the same result as 
we got with 12 variational parameters, and hence this 12 parameter form was used in the section below. 
This form gives dimensionless standard-deviation, defined as 



<<(5-5 a ) 2 )-«5-5 a » 2 



V <s 2 > 

being approximately equal to 21% . Here, 

J(X)d 3 R 1 d 3 R 2 d 3 R 3 



(X) 



J(l)d 3 R 1 d 3 R 2 d 3 R 3 



As this is not too small, in our main calculations we have made a minimal use of this further approx- 
imated area S a (only for the to-be-differentiated right y/J part (see eqUU] ) of the kinetic energy term 
and here only for derivatives of the exponent). 
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4 Solving the integral equations 



In eq.© for k = I = 1 (a diagonal term), we used the linear independence of y l7 Zj and R4 (see eq.©) 
to take xi(R-i) outside the integrations w.r.t. y 1 and z\. For the off-diagonal term with k — 1 and 
/ = 2 we replaced y 1 and Zi with R 2 and R 3 , with Jacobian of transformation as 8. For regulating the 
space derivatives of the exponent of / (see the three sentences immediately following eq. (TT5j) above) we 
temporarily replaced S in it by its quadratic approximation S a . As a result, we obtained the following 
equation: 



^ _ J-Vl,, + 24C a d 2 -\c-E c + 4mVi(Ri) 

R\ +Rl+2R\ 



J d 3 R 2 d 3 R 3 exp ^ - b s kfS^j exp 



2d 2 



7 gi exp 



6m(27rd 2 ) 3 



3(27Td 2 ) 3 



X2(R 2 ) =0, 



(17) 



with, written up to accuracy 4, 



gi = -1.4417 + 0.0258xi 2 + 0.0258a; 2 2 

0.0258y 2 2 + 0.0254y 3 2 + 0.0258zi 2 



0.0254a;3 2 



0.0258?/i 2 



0.0258z 2 2 + 0.0254z 3 2 . 



(18) 



For the consistency of £fc(y/c) and Cfe( z fe) with cq.© u) = 1/md 2 = 0.416GeV. For convenience in notation 
we take C = — C/3. Here in the first channel for k = 1 the constituent quark masses has been replaced 
by the reduced mass /Z12 = M1M2/ (Mi + M2), where Mi and M 2 are masses of hypothetical mesons; a 
similar replacement has been done in ref.[3]- 

At this stage we can ht C to a kind of "hadron spectroscopy" for our equal quark mass case: 
For the large separation there is no interaction between Mi and M 2 . So the total center of mass 
energy in the large separation limit will be the sum of kinetic energy of relative motion and masses of 
Mi and M 2 i.e. in the limit Ri — > 00 we have 



2/ii 



-V 2 



Mi + M 2 



Xi(Ri) = S c xi(Ri). 



(19) 



By comparing, in this limit, eq. (IT9l and eq. flTTl) we have Mi + M 2 = 4m + 3uj — 8(7/3. (A use of the first 
term of eq.© for the colour-basis diagonal matrix element of eq.© gives — 4C/3 = fiU! 2 /2 — fiui/2md 2 , 
giving 24C d 2 — 3w/2 for the reduced mass [i of a pair of equal mass quarks being m/2 ; the diagonal 
elements in any form of the / model for the gluonic basis are the same as those for the colour basis.) By 
choosing Mi + M 2 = 3w we have C = 3m/2 = 0.45GeV. This choice of the hypothetical meson masses is 
the one frequently used in ref. [3] for an illustration of the formalism; when we incorporate flavour and 
spin dependence [20] the same fit, something like in ref. [16], we plan to fit our quark masses to actual 
meson spectroscopy. We can then choose to fit even the parameter C or Co of our potential model to 
hadron spectroscopy rather than deciding it, as in ref. [3] and the present work, through a combination 
of baryon radii and harmonic oscillator model. But we do not see any reason why the qualitative effects 
(for example, an angle dependance, see the section below) pointed out through the present work should 
disappear for a phenomenologically explicit case. 

Completing our integral equations before finding a solution for two x's, for k — 2 in eq.© we took 
X2 (R-2 ) outside of integration for the diagonal term, for the off-diagonal term we replaced y 2 and z 2 by 
Ri and R3 and replaced S by S a . This resulted in 



1 

VI 



3a; 
2 



2// 



34 



24C a d 2 ~lc-E c + 4m )x 2 (R 2 ) 



d 3 i?id 3 i? 3 exp ( - b s kfS) exp 



Rj + R 2 + 2R 2 
2d 2 



6m(27rd 2 ) 3 



Si exp 



exp ( - -b s kfS a 



9(Wfr 4GK3 2C )~ 3(2^d 2 )3 



xi(Ri) = o. 



(20) 
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In the 2nd channel, for k = 2, the constituent quark masses are replaced by the reduced mass 
/X34 = M3M4/ '(M3 + M4), where M3 and M4 are masses of hypothetical mesons. 

Now we solve our two integral equations. As our space derivatives have been regularized, we no longer 
need further-approximated S a and we replace this by the original S in ea. (|20p . Below we take Fourier 
transform of eq. (|17|) . This gives us a nine dimensional integral of, amongst others, exp(—b s kfS). Non- 
separability of S did not allow us to formally solve the two integral equations for a non-trivial solution for 
Xi and \2 as in ref. [3] and we had to pre-specify a form for ^2(1^2) in eq. (|17| and of xi(Ri) in eq. (|20| . 
(As long as all the functions, including the meson wave functions and the gluonic field overlap factor /, 
are separate in Ri and R2 , we can everywhere replace Xi an( A X2 by their analytical integrals which 
themselves simply multiply if xi an d X2 do, can solve the resulting linear equations for these integrals 
and can write the T matrices and phase shifts directly in terms of these integrals. This is what is done in 
ref. [3] [16] , but it is hard to think how to generalize this very specialized technique to a case like us where 
the / factor is not separable in Ri and R-2-) Compelled to use, thus, Born approximation (something 
already in use [1] for meson-meson scattering; our numerical results mentioned below also justify its use 
here) for this we used the solutions of eqs. (fT7|) and (|20p in absence of interactions (say by letting kf 
approach to infinity, meaning / = 0) for xi(Ri) and %2(R2)- We chose the coefficient of these plane 
wave solutions so as to make Xi(Ri) a s Fourier transform of 8{P\ — P c (l))/P 2 (1) and X2(R2) as Fourier 
transform of S(P 2 - P c (2))/P 2 (2), with P c (l) and P c (2) defined below just after eq.([2S]). Thus we used 



X2(R2) = y^ ex P (*P2-R2) 



(21) 



inside the integral to get one equation (after a Fourier transform with respect to Pi and kernel e lPl - Rl ) 
as 



3w 



2/Xii 



-E C Xi (Pi) = 



Yf^jT J d 3 Rid 3 R 2 d 3 R 3 cxp|i(Pi.Ri +P 2 .R 2 )| 



exp ( — b s kfS^j exp ^ — 



R\ + R\ + 2Ri 
2d 5 

6m(27rd 2 ) 3yl ' 9(2^d 2 ) 3 V 3 J 3(2^d 2 ) 3 



32 / 
191 + n "2ird 2 Y V 



(22) 



with xi(Pi) being Fourier transform of xi(Ri)- The formal solution [3] of ea.([2^|) can be written as 



Xi(Pi) 



5(P - P c (l)) 



exp 



b s kfS ) exp ( — 



R\ + R\ + 2Rj 



2d 2 



^ + Y (-^Ri-2C)-"-( Ec -4m) 



(23) 



with 



Ai(Pi 



IL 

2/Xn 



3lu — E c — ie. 



If we choose s-axis along Pi and choose z-axis in such a way that ax-plane becomes the plane 
containing Pi and P2, the above equation becomes 



^(1) 



Ai(Pi)' 



(24) 



where, in the notation of eq.([T5 
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Fi — r, 5 ,« / dx\dx2dx 3 dy\dy2dy 3 dz\dz2dz 3 exp { iP(xi + X2 cos (9 + Z2 sin6l)l exp f — fogfefS*) 
167r 5 c( b y_ oc I ' J V / 

x\ + y{ + z\ + x\ + y\ + 4 + 2{x 2 3 + y 2 3 + z 2 ) 



exp 



2d 2 



6m 



.91 



+ y { - 4C(^ + y\ + z f ) - 2C\ --(E C - 4m) 



Here 9 is the angle between P2 and Pi and because of elastic scattering P\ = P2 = P. From eq. (j24]) we 
can write, as in ref. [3], the 1, 2 element of the T- matrix as 



T12 = 2/ii 2 -P c Fi. 



(26) 



Here P c = P c (2) = P c (l) = ■ S /2^{E C - {M x + M 2 )) and M x = M 2 = 3w/2; see paragraph after ea. ([T9")) . 
Using the relation 

s = I - 2iT = cxp(2iA) 



or 



i 


n 




• 


)-( 












P22 . 




v 





2i 



Su Sis 
S21 822 



between s matrix and the T matrix (actually in the form of elements Sij — —Tij for i,j = 1, 2) we got 
different results for phase shifts for different values of center of mass kinetic energy T c and the angle 9 
between P x and P 2 ; we have used the Born approximation to neglect higher powers in the exponential 
series. We also probed different values of the parameter kf. 

For a comparison, we also did the much less time consuming (but approximate) calculation using S a 
in place of S in eq. ([22|) . This allowed us separating the 9 variables dependence of the integrand as a 
product, resulting in three triple integrals to be only multiplied, making the convergence very fast in the 
numerical computation of the integral. Thus we had instead 



Xi(Pi) = 



Ai(Pi) 



F. 



with 



(27) 



F = 



1 



167T 5 d 6 



xi 2 + x\ + 2x\ 



J dx\ dx-2 dx-A J /i (xi , X2 , £3) exp 
b s kf(^ax\ + dx\X2 + JX1X3 + bx\ + gx2X 3 + cx^j + iP{x\ + X2 cos ( 
dyx dy 2 dy 3 < f 2 (yi , 2/2 , Vs, ) exp 



2d 2 



yi 2 +y 2 + 2y 2 



2d 2 



-b s kf (ay\ + e yi y 2 + ky x y 3 + by\ + hy 2 y 3 + cy\ 

/ dz 1 dz 2 dz 3 lf 3 (zi, z 2 ,z 3 ) exp 
J — 00 



>Q(x) x Q(z) 



zi 2 + z\ + 2zf 



2d 2 



-b s k f 



(az 2 + fz±z 2 + lziz 3 + bz\ + iz2Z 3 + cz 2 ^j + iPz 2 sin ( 



>Q(x) x Q(y) 



(28) 



Here 



fi(x u x 2 ,x 3 ) = - g|r(- 1-4417 + 0.0258a:i 2 + 0.0254X2 2 - 4.1914 x lO" 7 ^^ + 0.0258 
■ ( - 4Cxl - 2Cj - § (E c - 4m 



32 i 
9 
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Myi,V2, 2/3 ) = - ^ (0.0258 yi 2 + 0.0254y 2 2 - 4.1914 x lO" 7 ^ + 0.0258y 3 2 ) - kfCyl 
f 3 {z u z 2 , z 3 ) = - B | r (0.0258z 1 2 + 0.0254z 2 2 + 5.1396 x KT 6 ^ + 0.0258z 3 2 ) - ±f Cz\, 

Q(x) = dx\dx 2 dx 3 exp — 3:1 + 2 x Js +2a:3 — b s kf (ax\ + dx\x 2 + jx\x 3 + bx\ + gx 2 x 3 + cx^j + 
iP{x\ + X2 cos 9) , 

Q{y) = S-00 dy 1 dy 2 dy 3 exp - Vl + 2 ^ 2Vi - b s k f (ayj + ey Y y 2 + ky x y 3 + by\ + hy 2 y 3 + cyfj 

Q(z) = dz\dz 2 dz 3 exp — Zl + 2 ^ +2z3 — b s kf(azl + fziz 2 + lziz 3 + bz 2 +iz 2 z 3 + czfj +iPz 2 sin( 

For this choice of S, we also calculated the phase shifts that are reported in the next section. 
By treating eq. (|20p in the same fashion as that of eq. (|17[) and using the Born approximation 

Xi(Ri) = y|exp(iP 1 .R 1 ) (29) 

it was checked that the results for phase shifts remain same. Actually eq. (|20|) and eq. (fl~7| become identical 
if we interchange Ri and R2. 

5 Results and conclusion 

Fig[T]shows our results, with kt defined by eq. tfTTj) taken as 0.5, for the phase shifts for a selection of center 
of mass kinetic energies for different angles between Pi and P2 (Some numerical uncertainty appears at 

0. 15GeV for 9 = 0. When we further explored the region between 0.14GeV and 0.16GeV there appeared 
fluctuations in the results. For smoothness in graph we have neglected the data point at 0.15GeV in figfl] 
at 9 = and interpolated data points are taken there.) 

We found no numerical fluctuations for kinetic energies above 0.16GeV, and thus we conclude that 
in this kinematical range the scattering angle has large effect on phase shifts, indicating a true gluonic 
field effect. (The origin of this angle dependence is the exponent S which is essentially in a model of 
W12 Wilson loop, a pure gluonic theory expectation value; we do not get any angle dependence if this S 
is not used. So the angle dependence emerges from gluonic field related to the area law, and is thus a 
QCD effect.) By increasing the scattering angle the phase shifts become large. We noted that a faster 
convergence of the nine-dimensional integration (see eq. (|25[l ) for large kinetic energy values was possible 
for smaller values of the parameter kf, for a decrease in kf of 0.1 the CPU time reduced at least three 
times to that for the previous value. Thus we used the smaller value of kf = 0.5 mentioned in ref.[4] to 
get phase shifts for a larger set of kinetic energies resulting in smoother graphs. For the above mentioned 
value C = -0.0249281GeV 3 (meaning lo = 0.665707GeV and d = 0.441357 fm) used in ref.[]3] (giving the 
99.4% overlap of the wave functions) we found that, at T c = O.lGeV, there is about a 1 degree change 
in phase shift for a 30 degree change in scattering angle 9, even larger in magnitude than the phase 
shifts of figJT]for the corresponding T c given by our routine value d = 0.556 fm. So we can say that the 
characteristic feature of angle dependence will remain if we study the scattering of some realistic meson 
meson system by taking sizes of mesons accordingly and adjust the parameters of the Gaussian wave 
functions to simulate realistic linear plus Columbic potential eigenfunctions. 

For a comparison with the other crude forms of / previously used, we show in the following fig|5]the 
average of these kf = 0.5 phase shifts over our selection of angle 9 values along with the corresponding 
phase shifts for other forms of / i.e. exponent in / being proportional to S a , Ylj^-j r 1j an d zero; the phase 
shifts were found to be independent of the angle 9 for all these older forms of / and hence there was no 
need to take any angle-average for these other forms. This figure shows that in comparison to kf = 
(sum of two body potential model) we get relatively very small coupling with S, S a and Gaussian form 
in /. The introduction of a many-body interaction in the previous (Gaussian) form of / resulted in a 
reduced meson-meson interaction. In ref. [3] this reduction was noted as decreased meson-meson phase 
shifts. So there are less chances of making a bound state with modifications in sum of two body approach 

1. e. the inclusion of gluonic field effects significantly decrease coupling between two mesons in a q 2 q 2 
system. Phase shifts are much less than 1 radian which indicates the validity of Born approximation. 
The phase shifts we get are lesser than reported by others who have used Born approximation pQ but not 
used / factor in off-diagonal terms. 
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It is to be noted that the S a form also does not result in any angle dependence, although in contrast to 
the / = 1 and Gaussian form there is apparently no a priori reason to expect such an angle independence 
for the use of S a - This may be because S a is almost Gaussian with a little mixture of X1X3 and Z1Z3 
terms (see ea. ()16|i and the parameter values reported just below it) or because S a can be converted to 
a Gaussian form by a completing of squares. As for a comparison of the S a phase shifts with the angle- 
average of the S phase shifts, it can be pointed out that the height of phase shift with S a became less 
than that with original S but the shape remains identical. Perhaps this indicates that S a simulates well 
some variations resulting from the original S form. In figO if we compare graph of Gaussian form with 
the angle averaged phase shifts using S in / we find that as compared to Gaussian form the graph of 
other forms is closer to kf = 0, though the height of graph with kf — is still very large as compared to 
both Gaussian form of / and that of S in /; see figJT] which clarifies any ambiguity, if so, in figj2] about 
the kf=0 results. 

Fig El reports most of the results for the higher values 0.57 of kf mentioned in ref. [4 . This value 
is more precise for their form of model and the crude area expression in the exponent in it. But our 
numerical calculations for this value turned out to be more demanding and thus were done for a smaller 
selection of kinetic energies. The numerical uncertainties for this value are for 9 = -| and the kinetic 
energy between O.llGeV and 0.12GeV; the results for this value of 9 are in figEJ 

A value of kf = 1.0 higher than 0.5 and 0.57 mentioned above has been reported in ref. [21]. Although 
this work analyzes a relatively limited collection of geometries (only squares and tilted rectangles), we 
have tried to see effects of using a higher kf. The numerical problems for large kf implied in the above 
mentioned numerical convenience for smaller kf did not allow us to get results for kf=l in a manageable 
time even for T c = O.lGeV. The best we could do was to do a number of calculations for kf — 0.8; the 
resulting phase shifts from these are shown in fig[J]except for the phase shifts for 9 = ir/2 that are reported 
in figlHl showing some numerical uncertainties at the kinetic energy=0.13GeV and at 0.49 — O.SlGeV. Based 
on figH] and figlHl we expect that for higher values of kf the results will remain qualitatively same and 
do not expect any new feature to emerge for kf = 1. 

We did a partial wave analysis for the phase shifts in fig. [TJ For this we have projected our angle 
dependence on m = spherical harmonics. The angle dependence is independent of azimuthal angle <fi 
so the partial wave expansion will only contain terms independent of angle <f> or terms with m=0. This 
analysis shows that below 0.2GeV the S-wave is very much dominant as can be seen by a sharp rise of 
graphs towards lower side in fig. [5]and fig. linear 0.2GeV and it is also justified from fig. [T]which shows a 9 
independence below 0.2GeV i.e. purely an S-wave. Our partial wave analysis shows that only even partial 
waves are present. Furthermore partial wave phase shifts decrease as we go from S-wave to I-wave, as is 
clear from large reciprocals in fig. [9j and they become negligible as compared to corresponding S-wave 
phase shifts as we go beyond G-wave. So here in figs. (|8|9p only S/D and S/G ratios are plotted. The 
reason for the absence of odd partial waves is that our phase shifts are symmetric around 9 = \ and the 
product of an even and odd function is an odd function giving rise to a zero result after integration. Thus 
phase shift is different for different angle and partial wave analysis of this angle dependence indicates 
presence of I = 2, 4 spherical harmonics along with the angle independent 1 = 0. We have found this 
extra presence of D, G waves in angle dependence only when we use e - constant * area form of /. But when 
/ is spherically symmetric, i.e. old Gaussian form, then there appears no D, G waves. Thus D, G waves 
must be a property of e - constant * area _ So e -~ constant * area m ay couple an 1=0 meson- meson system to 
1=2, 4,... etc. systems. This means two possibilities 1) 1=0 meson-meson system may have t-matrix and 
phase shifts to 1=2, 4,... final state meson-meson systems and 2) 1=0 may couple to 1=2, 4,... states as 
intermediate states in a polarization potential [19], through e - constant * area _ 
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Tc(GeV) 

Figure 1: For kf — 0.5, the comparison of phase shifts for different values of 9 using S in /. The graph 
with points is for 9 = 0, with steps for 6 = ^ and with dots only is for 8 = j. The graph with lines is for 
9 = 5; here the data is equally spaced as in the other graphs but data points are joined. 




Tc(GeV) 

Figure 2: Comparison of different forms of /. Graph with lines only corresponds to kf = 0. (Here data 
points are O.OlGeV apart but they are joined for clarity. This graph is also shown separably in figtZJ 
Only four portions of the graph are shown here: upper peak corresponds to actual data minus 28 and 
lower peak corresponds to actual data plus 5.) Graph with dots corresponds to S a in /. Graph with steps 
corresponds to gaussian form of / for kf = 0.075 as defined by eq. (fl"2|) . Graph with points corresponds 
to average of phase shifts for different values of angle 9 for S in / with kf = 0.5 defined by ea. ([TT|) . 
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Figure 3: Comparison of phase shifts for different values of 8 using S in / with kf — 0.57. Choice of 
curve shapes is same as in fig. [TJ But the results for 9 — tt/2 are shown in fig. [5] and not here. 
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Figure 4: Comparison of phase shifts for different values of 9 using S in / with fc/ = 0.8. Choice of curve 
shapes is same as in fig. [3] The results for 9 — tt/2 are shown in fig. [6] below. 
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Figure 5: Phase shifts for kf — 0.57 at 9 — ^ using S in /. 




0.2 0.4 0.6 0.S 1 1.2 1.4 



Tc(Gev) 

Figure 6: Phase shifts for kf = 0.8 at 6 = § using S in /. 



16 



0.5 1 1.5 

Tc(GeV) 

Figure 7: Phase shifts for kf =0. 
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Figure 8: S/D ratios for different values of T c 
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Figure 9: S/G ratios for different values of T c . 
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